% Code for CC/MLI ->PC paired recordings batch processing. It is run from
% the analyze CC-PC_allPairs (this is for IPSCs in postsynaptic cell, for
% on cell on both check CC_MLI_Paired.m)

% Run this section if you need to average a bunch


fileName=' '

%% Use user input to determine which one is the presynaptic channel
formatpattern = fullfile(dir_to_search, ['*.ibw']);
dinfo=dir(formatpattern);

% determine which channel is presynaptic using variance of trace
tempDataStruc(1)=IBWread([dinfo(1).folder '\' dinfo(1).name]);
tempDataStruc(2)=IBWread([dinfo(1+numel(dinfo)/2).folder '\' dinfo(1+numel(dinfo)/2).name]);
variance1=var(tempDataStruc(1).y(1:20e3));
variance2=var(tempDataStruc(2).y(1:20e3));
if variance1>variance2
    preCh=str2num(tempDataStruc(2).bname(3));
else
    preCh=str2num(tempDataStruc(1).bname(3));
end

% define postCh
if preCh == 0
    postCh=1;
else
    postCh=0;
end


%OLD SECTION TO DO THE ABOVE MANUALLY, ANNOYING TO DO SO SO I REPLACED IT
%BUT KEEP THIS JUST IN CASE
% tempDataStruc(1)=IBWread([dinfo(1).folder '\' dinfo(1).name]);
% figure;
% plot(tempDataStruc(1).y);
% a=dinfo(1).name;
% 
% title(a(1:3));
% 
% prompt = 'select presynaptic channel (type 0 or 1)         ';
% preCh = input(prompt)
% close gcf
% 
% 
% if preCh == 0
%     postCh=1;
% else
%     postCh=0;
% end


clear dinfo formatpattern;
%%

formatpattern = fullfile(dir_to_search, ['**/ch' num2str(preCh) '_*.ibw']);
dInfoNS = dir(formatpattern);
[~, idx] =sort([dInfoNS.datenum]); % sort structure indices by date of creation, to avoid 10 coming before 2.

dinfo=dInfoNS(idx); % turn the not sorted structure to a sorted one.

nTrials=numel(dinfo)
%% initialize variables

spikes=[];  %initialize variable for adding presynaptic spike locations in each trial
IPSCAmplitudeBaselineAT=[]; % for all trials
  IPSCAmplitudeStimAT=[]; %for all trials
  isiIPSCBaselineAT=[]; %for all trials
  isiIPSCStimAT=[]; %for all trials


%% the loop
for trial=1:nTrials
    trial
    dataStruc(1)= IBWread([dinfo(trial).folder '\' dinfo(trial).name]);
    
    postName=dinfo(trial).name;
    postName(3)=num2str(postCh);   

    dataStruc(2)= IBWread([dinfo(trial).folder '\' postName]);
    
    Dx=dataStruc(1).dx;  %take any IBW to find sampling rate.
    Fs=1/Dx;
    
    
    rawTraces(:,1,trial)=dataStruc(1).y;
    rawTraces(:,2,trial)=dataStruc(2).y;
%         rawTraces(:,2,trial)=medfilt1(dataStruc(2).y,20);

    time=([1:numel(rawTraces(:,1))]*Dx)';
    
    %% remove step artifact from presynaptic channel
    
    buffer=100;
    artifact=[1e5 2e5];
    
    rawTraces(:,1,trial)=rawTraces(:,1,trial)-medfilt1(rawTraces(:,1,trial),100); %remove baseline from presynaptic channel, i.e. remove step but keep spikes.
    rawTraces([artifact(1)-buffer: artifact(1) + buffer],1,trial)= NaN;
    rawTraces([artifact(2)-buffer: artifact(2) + buffer],1,trial)= NaN;
    
  %  figure;plot(time,rawTraces(:,1,trial)); %fig to review spike trace
  %  after baseline substraction
%     %% plot the cumulative charge in pC
%     
%     % figure; plot(time,cumsum(rawTraces(:,2)))
%     cumCharge=cumtrapz(time,rawTraces(:,2,trial)-mode(rawTraces(:,2,trial)));
%     in=figure;
%     
%     subaxis(2,1,1);
%     plot(time,cumCharge,'linewidth',2);
%     xlabel('time (s)')
%     ylabel('cumulative charge (pC)') % it's cumulative integral
% %     title ([num2str(recNumb)])
%     title([dinfo(trial).name])
%     vline(artifact*Dx,'--r'); %label where the stimulus occurs
%     
%     
%     % Make PSTH
    [pks, locs]=findpeaks(rawTraces(:,1,trial)*-1,'minPeakprominence',30,'minpeakDistance',50); %finds the presynaptic spikes
%     subaxis(2,1,2);
%     plot(time,rawTraces(:,1,trial))
%     
% %     
% %     A=psth(locs,200,1/Dx,1,numel(rawTraces(:,1,trial)));
% %     xlabel('Time (ms)');
% %     ylabel('Average firing rate (spikes/s)');
% %     
    %% find instantaneous firing rate of presynaptic cell and plot in a separate
    % plot
%   
    A=locs(1:end-1);
    B=locs(2:end);
    isi([1:numel(B-A)])=B-A;
    hold on;
    iFiringRate=1./(isi*Dx);
%     figure
%     plot(B*Dx,iFiringRate([1:numel(B)]))
%     xlabel('time (s)')
%     ylabel('Instantaneous firing rate (spikes/s)')
%     vline(artifact*Dx)
% %     
    %% Have a moving window for calculating charge
    
    postBS=[];
    %%postBS=rawTraces(:,2,trial)-mode(rawTraces(:,2,trial)); % this sucks, use 5th percentile instead
    % create a baseline trace, using the 5th percentile in a window moving
     postBS=rawTraces(:,2,trial)-prctile(rawTraces(:,2,trial),5);
    
%     windMs=1000; %in miliseconds
%     windSam= windMs/(Dx*1e3);


%     for h=1:numel(rawTraces(:,2,trial))
%         if h<windSam
%             baselineTemp(h)=prctile(rawTraces([1:windSam],2,trial),5);
%         else
%             if h>numel(rawTraces(:,2,trial))
%               baselineTemp(h)=prctile(rawTraces([end:end-windSam],2,trial),5);
%         else  
%             baselineTemp(h)=prctile(rawTraces([h-round(windSam/2):h+round(windSam/2)],2,trial),5);
%             end
%         end
%     end
%     postBS=rawTraces(:,2,trial)-prctile(rawTraces(:,2,trial),5);


% % this is very very slow and works OK but removes any DC effects (tonic?)
% postBaseline = RankOrderFilter(rawTraces(:,2,trial), windSam, 2);
% postBS=rawTraces(:,2,trial)-postBaseline


%OLD
%     windMs=1000; %in miliseconds
%     windSam= windMs/(Dx*1e3);
%     for i=1:(numel(postBS)-windSam); % remove window pad on both sides
%         chargeRate(i)=trapz(postBS(i:i+windSam))*Dx;
%     end
%     
%     time2=time(1+windSam/2:end-windSam/2);
%     figure;
%     plot(time2,chargeRate(:));
%     
    
    %% Bin trace for real binned charge without moving window (avoid delay)
    
    binWidth2=0.5; %in seconds
    binWidth2Sam=binWidth2/Dx;
    for i=1:(numel(postBS)-1)/binWidth2Sam
        chargeRate3(i)=trapz(postBS(1+(i-1)*binWidth2Sam : 1+binWidth2Sam+(i-1)*binWidth2Sam))*Dx;
        time3(i)=time(1+(i-1)*binWidth2Sam)+binWidth2/2;
    end
  
    %% find the firing rate for each time bin.
    
    bounds=[5.5 9.5]; %start and end points, 5.5 is 5(start) + windowMs/2 (1s/2))
    timeBin=0.5; % window for calculating firing rate and charge rate.
    spikeTimes=B*Dx; %convert from indices to seconds.
    
    %find the average firing rate and charge rate for each time bin
    nBins=((bounds(2)-bounds(1))/timeBin )+1;
    for i=1:nBins
        midPoint=bounds(1)+(i-1)*timeBin;
        startPoint=midPoint-timeBin;
        endPoint=midPoint+timeBin;
        
        firstSpikeIdx=find(abs(startPoint-spikeTimes)==(min(abs(startPoint-spikeTimes))));
        lastSpikeIdx=find(abs(endPoint-spikeTimes)==(min(abs(endPoint-spikeTimes))));
        meanFR(i)=mean(iFiringRate(firstSpikeIdx:lastSpikeIdx));
        
        % meanChargeRate(i)=mean(chargeRate(((startPoint-0.5)/Dx):((endPoint-0.5)/Dx))); %Substract the timeBin to correct for padded window for calculating chargeRate (i.e. first 500 and last 500 ms are gone)
        
    end
    
   % baselineChargeRate=mean(chargeRate(1:round(4.5/Dx)));
   % percChargeRateReduction(:)=1-(meanChargeRate(:)./baselineChargeRate);
    
%     figure; plot(meanFR(:),percChargeRateReduction(:),'k.','markersize',30)
%     set(gcf,'color','white');
%     set(gca,'FontSize', 18);
%     box off

%% IPSC Analysis section
%find IPSCs (size and location)
[pksIPSC, locsIPSC]=findpeaks(postBS,'minPeakprominence',20,'minpeakDistance',50);

% figure; %figure for validation of IPSC detection
% plot(time,postBS)
% hold on
% plot(locsIPSC*Dx,pksIPSC,'r.','markersize',10)
% for i=1:numel(locsIPSC)
%     plot([locsIPSC(i)*Dx locsIPSC(i)*Dx],[0 pksIPSC(i)]);
% end

%sort into baseline and stim times
locsBaselineIPSC=locsIPSC(locsIPSC<artifact(1)); %find ipscs that ocurr before stim onset
pksBaselineIPSC=pksIPSC(locsIPSC<artifact(1));

locsIPSCStim=locsIPSC(locsIPSC>artifact(1) & locsIPSC<artifact(2)); %ipsc that ocurr between the two artifacts, i.e. during stimulation time
pksIPSCStim=pksIPSC(locsIPSC>artifact(1) & locsIPSC<artifact(2)); %same for peask

%convert IPSC locations into intervals in miliseconds
isiStimIPSC=(locsIPSCStim(2:end)-locsIPSCStim(1:end-1))*Dx*1e3;
isiBaselineIPSC=(locsBaselineIPSC(2:end)-locsBaselineIPSC(1:end-1))*Dx*1e3;


  %% Append data for each loop iteration
  %presynaptic spikes
  spikes=[spikes ; locs]; % stack all spikes together for psth
%  chargeRateAllTrials(:,trial)=chargeRate;
  BinnedFRAllTrials(:,trial)=meanFR;
  %BinnedChargeRateAllTrials(:,trial)=meanChargeRate;
  %fractionalChargeRateReduction(:,trial)=percChargeRateReduction;
%   postSynaptic IPSC
  IPSCAmplitudeBaselineAT=[IPSCAmplitudeBaselineAT ; pksBaselineIPSC];
  IPSCAmplitudeStimAT=[IPSCAmplitudeStimAT ; pksIPSCStim];
  isiIPSCBaselineAT=[isiIPSCBaselineAT ; isiBaselineIPSC]; %names are confusing but i'm just appending each trial to a general variable
  isiIPSCStimAT=[isiIPSCStimAT ; isiStimIPSC]; % same as above.
  realBinnedChargeRate(:,trial)=chargeRate3;
end
    

%% Make PSTH of all trials (presynaptic firing)
 A=psth(spikes,200,1/Dx,nTrials,numel(rawTraces(:,1,trial)));
 FR_Temp=A.Vertices(:,2);
 FR_Temp=FR_Temp(3:end); %trim ends
 FR_Temp=FR_Temp(1:end-4); %same as above
 FR_Temp=FR_Temp(1:5:end);
% xlabel('Time (ms)');
% ylabel('MLI firing rate (spikes/s)');
% box off
% set(gcf,'color','white');
% set(gca,'FontSize', 18);

%plot IPSC charge binned for all trials.
% figure;
% hold on
% for i=1:nTrials
%    
% plot(time3,realBinnedChargeRate(:,i),'color',[0 0 0 1/4])
% end
% plot(time3,mean(realBinnedChargeRate,2),'color',rgb('firebrick'),'linewidth',2)
% Oylim=ylim;
% ylim([0 Oylim(2)])
% set(gca,'FontSize', 18);
% set(gcf,'color','white');
% box off
% xlabel('time (s)')
% ylabel('IPSC Charge (A.U)')
% vline([5 10],'--k');

meanChargeBinned=mean(realBinnedChargeRate,2);
baselineChargeMean=mean(meanChargeBinned(1:10));
ccFiringChargeMean=mean(meanChargeBinned(11:20));
effectSize=(ccFiringChargeMean-baselineChargeMean)/baselineChargeMean 


fclose all


